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Abstract - We investigate the stability of synchronization in networks of delay-coupled excitable 
neural oscillators. On the basis of the master stability function formalism, we demonstrate that 
synchronization is always stable for excitatory coupling independently of the delay and coupling 
strength. Superimposing inhibitory links randomly on top of a regular ring of excitatory coupling, 
which yields a small-world-like network topology, we find a phase transition to desynchronization 
as the probability of inhibitory links exceeds a critical value. We explore the scaling of the critical 
value in dependence on network properties. Compared to random networks, we find that small- 
world topologies are more susceptible to desynchronization via inhibition. 
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Introduction. — Studies of complex networks have 
iparked tremendous scientific activities in many research 
areas and the analysis of network topologies in real- world 
^stems has become a field of large interest. For instance, 
there is evidence that neuronal networks on the level of sin- 
gle neurons coupled through synapses or gap junctions, as 
well as on the level of cortex areas and their pathways ex- 
hibit the small- world (SW) properties [TP] . The high clus- 
tering coefficient of the SW networks enhances local com- 
ipunication efficiency, while the small shortest path length 
enables efficient global communication [3]. Thus, the SW 
Architecture is optimal for processing and transmission of 
signals within and between brain areas. However, the syn- 
qhroniz ability of small-world networks depends in a deli- 
cate way upon the network topology [3]. Next to this 
structural aspect, inhibition plays a prominent role in 
many neural processes [S]. Without an inhibitory mecha- 
nism, excitation in a compound system would not decay, 
but spread through the whole network, finally leading to 
persistent spiking of all neurons. Thus, encoding and pro- 
cessing of information would be impossible. 

In this Letter, we combine both fundamental aspects 
- inhibition and SW property - in order to emphasize 
the important interplay of excitation and inhibition in 
complex networks. We start with a regular ring network 



that consists of purely excitatory links with delay. Thus, 
it exhibits strong and stable synchronization. Depend- 
ing on the initial conditions, both isochronous and clus- 
ter synchronization are possible implying multistability. 
Additional inhibitory connections, which we include in a 
SW-like manner [Uini , result in a loss of synchronization. 
A similar transition was reported for phase oscillators in 
Ref. |7j for unidirectional rings, but the effect of inhibi- 
tion upon excitable systems could not be treated by that 
model. For the node dynamics, we consider a generic 
model to demonstrate the fundamental relevance and im- 
portance of our findings in the field of neuroscience. In this 
area, synchronization can be related to cognitive capaci- 
ties j8] as well as to pathological conditions, e.g., epilepsy 
[HI . A better understanding of the loss of synchronization 
will eventually lead to future therapeutic treatments |10) . 
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Throughout this Letter, we consider a network of N 
delay-coupled FitzHugh-Nagumo (FHN) oscillators. The 
FHN system describes neuronal dynamics by a two- 
variable model [TT]. Because of its simplicity it can be 
considered as a paradigmatic model of excitable systems, 
which also occur in several other natural contexts ranging 
from cardiovascular tissues to the climate system P^I13| . 
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Fig. 1: (Color online) Dynamics in the synchronization man- 
ifold in dependence on the coupling strength C and delay r. 
The gray scale (color code) indicates the period of spiking os- 
cillations T, the black region corresponds to fixed-point dy- 
namics. Left and right insets show time series of the activator 
Us for (C — 0.3, T = 1) and (C — 1.5, r = 1.8), respectively. 
Parameters: a — 1.3, e = 0.01. 



Here, the network dynamics is described by 



N 



(lb) 



u, == u, + a, 



where Ui and Vi denote the activator and inhibitor vari- 
ables of the nodes i = 1, . . . ,N, respectively. The param- 
eter a determines the threshold of excitability. A single 
FHN oscillator is excitable for a > 1 and exhibits self- 
sustained periodic firing beyond the Hopf bifurcation at 
a = 1. Here, we will focus on the excitable regime with 
a — 1.3. The time-scale parameter e is chosen as e = 0.01. 
C is the coupling strength. G = {Gij}, i,j = 1,...,A^, 
denotes the coupling matrix that determines the topol- 
ogy of the network. In the following we will assume unity 
row sum of G. This ensures that each neuron receives the 
same input if the network is synchronized. The delay time 
T takes into account the finite propagation speed of an ac- 
tion potential. We investigate complete synchronization 
with (w,(i),w,(i)) = (w,(i),w,(i)) =x,(t) fori = 1,...,^, 
which is also known as zero-lag or isochronous synchro- 
nization. This state is a solution of qs. ([Ij and reduces 
the system's dynamics to 



x, = F(x3) + CU [x,(t - r) - x,(t)] 



(2) 



and the matrix H = {^l,^ ?.] 



withF(x) = (lu-u'/3-v]/e\ 
^ ^ \ u+a j 

The 2(A^ — 1) constraints of complete synchronization de- 
fine a two-dimensional synchronization manifold (SM) in 
the 2A^-dimensional phase space. 

As we operate in the excitable regime, the dynamics on 
this SM, in particular the period, will depend on the choice 
of the coupling parameters C and t as depicted in fig. [IJ 



The grayscale (color code) corresponds to the period T of 
the oscillations on the SM, which we find to follow T = 
r+^ with (5 ^ T accounting for a short activation time |14| . 
For small coupling strength C the incoming signal is not 
sufficient to trigger oscillations (black region). For small 
delay times r consecutive spikes run into the refractory 
phase of the previous one, which prevents oscillations as 
well. From here on we consider C and r sufficiently large 
such that the coupling induces oscillations. 

Stability analysis. — In the following we address the 
question whether the oscillatory solution on this mani- 
fold is transversely stable. The master stability function 
(MSF) |15| allows us to quantify this transversal stability. 
It can be calculated as largest Lyapunov exponent A from 
eq. (HI) linearized around eq. ([2]): 

C(i) = [i?F(x,(t)) - CH] C(t) + (a + */3)HC(t - r). (3) 

Here, £>F denotes the Jacobian of F. The idea of the MSF 
is to calculate the stability of a synchronized solution for 
an arbitrary topology matrix G. For this purpose, the 
parameter a-^iji represents a continuous parametrization 
of {Cvi\, where Vi, i = 1, . . . ,iV, are the eigenvalues of 
G. In the same sense, the vector ^ is a generalization of 
the variational vectors transformed to the corresponding 
eigensystem. In the (a,/3)-plane the MSF typically gives 
rise to regions with negative A. If all rescaled transversal 
eigenvalues {Ci^i} of a given network are located within 
this stable region, perturbations from the SM will decay 
exponentially and the synchronized dynamics will be sta- 
ble. Due to the unity row sum condition, G will always 
have one eigenvalue i^i = I. This longitudinal eigenvalue 
is associated with perturbations within the SM and is 
not relevant for the stability of synchronization. A(Ci'i) 
determines the type of dynamics on the SM. For peri- 
odic dynamics in the SM, as in the present case, we have 
A{Cui) = 0. 

Figured] depicts the MSF for the network of FHN oscil- 
lators given by eqs. ([T|). Dark (blue) colors mark the stable 
region. As an illustration the rescaled eigenvalues {Cfi} 
of a bidirectionally coupled ring (TV = 8) are shown as red 
symbols. The corresponding coupling matrix is given by 
Gi,i+i mod N = Gi^i-i mod jv = 1 (« = 1, • • ■ , A^) and zero 
otherwise. The rescaled longitudinal eigenvalue Ci^i = C 
is depicted by a black (red) square. All rescaled transver- 
sal eigenvalues (black (red) circles) lie inside the stable 
region indicating that the synchronization of the bidirec- 
tionally coupled ring is stable. 

Shape of stability region. — The MSF must be cal- 
culated for each combination of C and r. Although dif- 
ferent C and r lead to quantitatively different Lyapunov 
exponents A, the shape of the stable regions remains qual- 
itatively very similar. In particular, it is in very good ap- 
proximation given by the circle 5'((0, 0), C) with center at 
the origin and radius C (dotted circle in fig. [5]) indepen- 
dent of the specific values of C and r. The rotational sym- 
metry has recently been proved generally for large r |16j . 
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Fig. 2; (Color online) (a) Master stability function for a 
network of FHN systems given by eqs. ([T}. Dotted curve: 
S'((0,0),C). Red circles (square): Rescaled transversal (longi- 
tudinal) eigenvalues Cui of a bidirectionally coupled ring with 
N = 8 nodes. Parameters: a = 1.3, e — 0.01, C = 0.3, r = 1. 
(b) Scheme of a bidirectional regular network (A*' = 20, fc = 2), 
and (c) a random network (A^ = 20, fixed number of links 
kN) with excitatory coupling (gray (green) arrows) on which 
inhibitory links (black (red) arrows) are superimposed. 



Only for small t and C the stable region is slightly larger 
than the circle and shows a bulge around a = — C, /3 = 
|17| . The positive a-axis is always intersected at a = C, 
which corresponds to A(Ci'i) = as discussed above. For 
any choice of r and C that leads to periodic dynamics on 
the SM, the circle S{{0, 0), C) serves as a lower bound for 
the stability boundary. See the appendix for an analytic 
derivation of this circle S'((0,0),C) in the limit of large 
coupling strength and as a lower bound for all coupling 
strengths. We conclude that the stability of the synchro- 
nized periodic dynamics, if such a solution exists, depends 
only on the topology and neither on the coupling strength 
nor on the delay time. 

Excitatory coupling. — For excitatory coupling, i.e., 
Gij > 0, all eigenvalues of G are located inside the sta- 
ble region. Using Gershgorin's circle theorem [TB], which 
gives an upper bound of the eigenvalues, and the constant 
row sum assumption, all Gershgorin circles (i = 1, . . . , N), 
centered at Gu with radius X^iV^j ^^ ~ ^ ~ ^" because 
of the unity row sunr, lie inside the unit circle. Thus all 
rescaled eigenvalues {Cui}, i = 1, . . . ,N, are located in- 
side S'((0,0),C), i.e., inside the stable region. Networks 
with purely excitatory coupling will always exhibit stable 
synchronization . 

Inhibitory coupling. — As a consequence of this 
result, desynchronization can only be achieved by intro- 
ducing negative entries in the coupling matrix G, i.e., 
inhibitory coupling between neurons. This inhibition is 
a crucial feature in neural processes, e.g., to overcome 
unwanted synchronization associated with pathological 
states. 

Particularly, we consider the following variation of the 
Watts-Strogatz SW network [Tlin|: (i) Start with a one- 
dimensional ring of TV nodes, where every node is con- 
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Fig. 3: (Color online) Fraction of desynchronized networks / 
vs the probability of additional inhibitory links p for A'^ = 100. 
k varies from 6 to 30. Thin black curve: Example fit to /(p) 
(pc = 0.20387, h = 186) for k = 24. Number of realizations: 
500 for each value of k. Parameters as in fig. \T\ 



nected by excitatory links to its k neighbors on either 
side, (ii) For each of the kN links of the network add 
an inhibitory link with probability p connecting two ran- 
domly chosen nodes, (iii) Do not allow self-coupling or 
more than one link between any pair of nodes, (iv) Nor- 
malize the entries of the coupling nratrix G by dividing 
each row by the absolute value of its row sum. In the case 
that the row sum of the ith row is negative we set Gu = 2 
to ensure unity row sum. Figure[21[b) illustrates such a SW 
network for A^ = 20 and fc = 2, where gray (green) and 
black (red) arrows indicate excitatory and inhibitory cou- 
pling, respectively. For each realization of such a network, 
we determine the stability of synchronization by checking 
whether the full eigenspectrum Ci^i of the coupling ma- 
trix is contained in the stable region ^((O, 0), C). Hereby 
we compute the fraction / of desynchronized networks. 
Figure [3] shows / as a function of p for different coupling 
ranges k. This Figure is virtually identical for all delay 
times, that is, for all parameters within the color shaded 
area of Fig.l. To obtain this Figure we made use of the 
circular shape of the stable region of the master stability 
function. Only for very small delays or coupling strength, 
the stable region is slightly larger and thus the shape of 
the curves shown in Figure [3] might be shifted slightly to 
larger values of p. 

For fixed k a steep transition between synchronization 
and desynchronization takes place as p approaches a crit- 
ical value pc- This critical value Pc and the steepness 
b of the transition can be fitted with a sigmoidal func- 
tion f{p) = l/[e-'''^P-P'--'> + 1]. Figure H^a) depicts the 
critical probability pc for f{pc) = 0.5 in dependence on 
k/N for different network sizes. It can be seen that 
for SW networks pc follows a linear relation pc{k/N) = 
1.16fc/A^ — 0.07 independently of the network size N . Fig- 
ure UJ^c) shows the steepness 6 as a function of the network 
size demonstrating that the transition becomes increas- 
ingly sharp as N increases. This indicates a first-order 
nonequililibrium phase transition in the thermodynamic 
limit [231 . 

To verify whether this phase transition and especially its 
independence of the network size is common in networks 
with inhibitory links or unique to the SW structure, we 
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Fig. 4: (Color online) Critical value pc for different network 
sizes, (a) in dependence on k/N, (b) in dependence on k: SW 
networks A'' = 60 (dark (blue) circles), A'' = 100 (dark (purple) 
squares), and A'' = 500 (turquoise triangles). Random networks 
TV = 60 (black (red) crosses), A*' — 100 (light (orange) circles), 
and A'^ = 500 (lightgray (yellow) squares), (c) Steepness b for 
SW (circles) and random (squares) networks vs. N for constant 
k/N — 0.1. Inset: Blow-up for random networks, (d) pc vs. A'^ 
for A; = 10 for a SW (black (red) filled circles) and a random 
(black (red) empty squares) network. Number of inhibitory 
links (I) vs. A'' for constant k for a SW (gray (green) empty 
circles) and a random (gray (green) filled squares) network. 
Number of realizations: 500. 



construct a different network for comparison: The regu- 
lar excitatory network is replaced by a random network 
with fixed number of excitatory links kN equivalently to 
the regular network used before, as shown in fig. [D^c). 
We only consider realizations where this underlying ex- 
citatory network is fully connected. The construction of 
the inhibitory links then follows steps (ii-iv) as above. We 
find that a phase transition to desynchronization still oc- 
curs with critical probabilities of inhibitory links pc as de- 
picted in fig. m^a) by black (red) crosses, gray (orange) 
circles, and lightgray (yellow) squares for N = 60, 100, 
and 500, respectively. We observe, however, that the val- 
ues of Pc are higher, i.e., the random network can tol- 
erate more inhibitory links than the SW network before 
desynchronizing. Furthermore, the iunction pc{k/N) is no 
longer independent of the network size A'^. Instead, Pc is 
a function of log(A:) as can be seen in fig. Hljb), where 
Pc is plotted in dependence on k for the different net- 
work sizes N. Figure |DJd) depicts Pc in dependence on 
N for constant k = 10 for a random (black (red) empty 
squares) and for a SW network (red (gray) circles). For 
random networks, pc is independent of N for sufficiently 
large A^, while for SW networks it approaches zero. Re- 
call that Pc is the mean value of the ratio of inhibitory to 
excitatory links. Thus, we conclude that in SW networks 
with increasing network size TV but same local structure 
(constant k) an infinitesimally small ratio of inhibition to 
excitation is needed to prevent synchronization, while in 



a random network even for very large networks only a 
non- vanishing ratio impairs synchronization. We find in a 
SW network with constant k that the mean value of the 
number of inhibitory links (/) := PckN causing desyn- 
chronization scales as (/) = 1.16fc^ — 0.07A:A^ for small 
N and approaches zero for large A^ (see fig. IDJd) green 
(gray) empty circles). In contrast, in a random network 
(/) is proportional to A^ (see fig. Hljd) green (gray) filled 
squares), i.e., an increasing number of inhibitory links is 
needed. This difference to SW networks can be under- 
stood in an intuitive way: In a SW, any added inhibitory 
link is part of a shortest path for many pairs of nodes, 
as it shortens the mean path length considerably with re- 
spect to the underlying regular ring. In the random net- 
work, however, where the mean path length is relatively 
low even without added shortcuts, only few node pairs will 
gain shorter paths by adding inhibitory links. Consider- 
ing the dynamics on a network, perturbations from the 
synchronized state spread along the shortest paths first, 
changing the response of the receiving node, and informa- 
tion flow along longer paths will reach the receiving node 
only at a later time and will not influence the change of 
the initial response. In conclusion, if a large fraction of 
the inhibitory links is part of the shortest paths - like in 
the small-world topology superimposed to a regular ring - 
these inhibitory shortcuts become dominant. 

Conclusions. — We have shown how the interplay of 
excitatory and inhibitory couplings leads to desynchro- 
nization in networks of neural oscillators. The desynchro- 
nization is achieved via a phase transition from a com- 
pletely synchronized state. This can be seen as a first step 
towards an understanding of the robustness of different 
states of synchrony, e.g., cluster synchronization, in arbi- 
trary networks with weighted links or distributed delays. 
Note that for appropriate network topologies the frame- 
work of the MSF presented above can indeed be extended 
to cluster synchronization where the oscillators synchro- 
nize in M clusters with a constant phase lag 2tt/M be- 
tween subsequent clusters [12]. The corresponding SM is 
2M dimensional. Hence, M longitudinal eigenvalues ex- 
ist. The MSF, however, is again very well approximated 
by the circle S'((0, 0), C) and thus, we observe multistabil- 
ity between zero-lag and cluster synchronization. 

Excitable systems can be classified into type-I and type- 
II excitability [iniHn]- In addition to the generic type-II 
FitzHugh-Nagumo model used in this paper, we have con- 
sidered the normal form of a saddle-node bifurcation on 
an invariant circle (SNIC) as a generic model of type-I 
excitability [H]. For sufficiently large delay times and 
coupling strength the MSF is again given by the cir- 
cle ^((0, 0),C) implying that the previously obtained re- 
sults persist. In particular, the same phase transition oc- 
curs. This indicates that the phenomena observed here 
are generic for any excitable system. 

Appendix: Analytic approximation of the stabil- 
ity region. — The numerical calculation of the master 



p-4 



Loss of synchronization in complex neuronal networks with delay 



stability function has shown that S'((0, 0),C) is a lower 
bound and a very good approximation of the stable region 
for all T and C. As t and C increase the approxima- 
tion becomes even better. A Taylor expansion as done 
in Ref. [35] for the investigation of time-delayed feedback 
control of an unstable periodic orbit gives analytic insight 
in the problem. This analysis is very general and does 
not use the specific form of the local dynamics in terms 
of the FHN model. It only assumes that the synchro- 
nized dynamics is oscillatory with period T . Using a Flo- 



where ab and /3t, denote the values of a and /3, respectively, 
on the bounder of stability. Finally we obtain 
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-Pi 



C2 



a ■ 



X 



(11) 



Obviously 



/3i 



that 
For 



quet ansatz 
Q(i) = Q(i 



c = 



^ p{A+in)t 



Q(f) with the periodic function 



T) in Eq. © yields 

(A + ir!)Q(f) + Q(t) (4) 

= {D¥ - CH)Q(i) + (a + z/3)e-(^+'")^HQ(t - r). 

K + i^ is the Floquet exponent, whose real part coincides 
with the Lyapunov exponent in the case of a periodic orbit. 
Assume T = t. In the case of the FHN system this is an 
approximation since the period of the oscillations differs 
by a small activation time S <^ t from the delay time r 
following T = T + S. Then Q(f — r) can be substituted by 
Q(i): 

{A + in)Q{t) + Q{t) (5) 

= pF)Q(t) + [-C +{a + ^/3)e-(^+*^^)^ HQ(i). 



We expand the solution r(K) = A + ifJ of the eigenvalue 
problem defined by Eq. ([S]) in a Taylor approximation: 



Tin) 



0{k'). 



r(o) + r'(o)K 

Using r(0) = X + iuj and r'(0) = x' + ix" we obtain 



(6) 



(7) 



Note that k = if (a, /3) = (C, 0) corresponding to the 
dynamics within the synchronization manifold. Thus the 
first term in the Taylor approximation corresponds to the 
Goldstone mode, i.e., X + icu = for (a,/3) = (C, 0). Equa- 
tion ([7]) then becomes 



A + in^x'[-C+ia + il3)e 



-(A+in)T 



Here, we assume x" = 0. Separating Eq. 
imaginary part leaves us with 



]• (8) 

into real and 
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n 



x'{- 

x'e 



-C 

At I 



-At 



a cos(nT) + j3 sin(f]r)]}. 



-a sin(J7r) -f /3 cos(f7r) 



(9) 



Equation ^ can be solved numerically yielding the cir- 
cular stability region. On the border of the stability the 
real part of the Floquet exponent vanishes. Using A = 
in Eq. ^ yields after algebraic manipulations: 



Oib = 



— r2sin(r2r) 
X' 

Ocos(r2T) 

X' 



+ Ccos(r2r), 
Csin(fiT), 



(10) 



> C^ holds, demonstrating 
S'((0,0),C) is a lower bound for the stable region, 
large C the term C^ on the right hand side dominates. 
Thus, the boundary of stability is very well approximated 
by S'((0,0),C) for large coupling strength. 
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